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We propose an analytic response theory for the density matrix renormahsation 
group whereby response properties correspond to analytic derivatives of density ma- 
trix renormahsation group observables with respect to the applied perturbations. 
Both static and frequency-dependent response theories are formulated and imple- 
mented. We evaluate our pilot implementation by calculating static and frequency- 
dependent polarisabilities of short oligo-di- acetylenes. The analytic response theory 
is competitive with dynamical density matrix renormahsation group methods and 
yields significantly improved accuracies when using a small number of density matrix 
renormahsation group states. Strengths and weaknesses of the analytic approach are 
discussed. 



I. INTRODUCTION 



The density matrix renormalisation ffroup method 111 is now established as a powerful 
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tool for "difficult" electronic structure problems in physics and chemistry [2, 13|, 14, 15[. In 
molecular systems, it has been used to describe multireference correlation in medium-sized 
active spaces (20-30 active orbitals) for small molecules with complex bonding 
well as a local multireference correlation method in extended long-chain molecules, e.g. to 
describe excited states in conjugated molecules, using large active spaces of up to 100 active 



orbitals 



io|. 



Response properties, which represent the change in an observable as a function of an 
applied perturbation, are of interest in many physical and chemical applications. For ex- 
ample, geometry optimisation and vibrational frequencies both require the response of the 
energy with respect to changes in the nuclear coordinates, quantities usually known as nu- 
clear derivatives. Nuclear derivatives are examples of static response properties because the 
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perturbation does not depend on time. It is also common to consider frequency- dependent 
(i.e dynamical) response properties where the applied perturbation is a function of time. 
The most common time- dependent perturbations are fluctuating electric and/or magnetic 
fields. In extended systems, the frequency dependence of the response gives insight into the 
elementary excitations of the system and this can be used characterise the nature of the 
electronic ground-state 

In many electronic structure methods, response properties are obtained by so-called "an- 
alytic" techniques. Analytic response theories of this kind at linear and higher orders have 
been developed and implemented for most electronic structure methods, including Hartree- 
7ock [12], density functional 13], coupled cluster multi-configurational self-consistent 
isj l, and Moller-Plesset perturbation theories 
of these developments may be found in Ref. 



161] . A review of the formal theory and some 



17|. The name "analytic" is used because the 



response properties evaluated (e.g. the perturbed energies) correspond strictly to derivatives 
of the ground-state energies or quasi-energies [l3, Q, Q evaluated in the presence of the 
perturbation, using the same level of approximation for the (quasi-) energy with and without 
the perturbation. 

In contrast, response properties in the density matrix renormalisation group have typically 
been obtained using a quite different approach that appears natural within the DMRG. In 
the DMRG, the wavefunction is expanded in a set of many-electron states that are adapted 
to the state of interest. To obtain a response property, one can choose to solve response 
equations using basis states that are adapted not only to the zeroth order state but also to 
the calculation of the state's response. These response methods, which have proven very 



useful in the calculation of dynamical response in DMRG model Hami 
by the name of Lanczos- vector DMRG [20] , correction- vector DMRG 
DMRG 



tonian calculations, go 
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22l | , and dynamical 



23j. More recently, explicit real-time propagation o 



the DMRG equations has also 



been used to obtain high-frequency response properties 2J]. A recent review of all these 
DMRG response methods can be found in Ref. 25 1. 

In the current work we return to an analytic formulation of response theory within the 
density matrix renormalisation group, in a way that parallels the description of response 
properties in other electronic structure methods. We use as our starting point the wave- 



function based (matrix-product state) formulation of the DMRG 



m 



27|. As we shall 



see, the analytic response approach has a number of strengths and weaknesses compared 
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FIG. 1: One-site DMRG block configuration. L" tensors are associated with the left block, R" 
tensors with the right block, and the middle site is site p in Eq. [TJ 

to earlier DMRG response methods. To understand these strengths and weaknesses better, 
we perform a series of benchmark static and frequency- dependent polarisabiUty calculations 
on oligo-diacetylenes that compare the behaviour of the earlier dynamical DMRG method 
with our analytic response DMRG approach. Using our data we examine the scaling of the 
polarisability as a function of the number of monomer units. 

II. TIME-INDEPENDENT AND TIME-DEPENDENT DENSITY MATRIX 
RENORMALIZATION GROUP EQUATIONS 

The density matrix renormalisation group works with a variational ansatz for the wave- 
: unction The simplest ansatz to analyse is the "one-site" form of the DMRG wavefunction 



271, 



281]. For the block- configuration depicted in Fig. [H the wavefunction takes the form 



1^) = ^L"i...L"''-iC"''R"''+i...R"*|ni...nfe) (1) 

{n} 

The L" and R" renormalisation tensors satisfy the orthogonality conditions 

J2 L'^^L" = 1 (2) 

n 

J2 = 1 (3) 

n 

and formally define the sequence of renormalisation transformations to obtain basis states 
{/}, {r} for the left and right blocks in Fig. [H (Note that in Eqs. ([2]), ([3]) we have 
dropped the sub-indices on n as these conditions are not specific to any given site. We will 
use a similar convention throughout to avoid a proliferation of unnecessary indices). The 
coefficient tensor C" gives the expansion coefficients of the wavefunction in the superblock 
basis {/} (S> {np} ® {r}. When viewed as a flattened vector c it satisfies the normalisation 
condition c^c = 1. 

The DMRG energy is minimised when the tensors satisfy certain equations. For the 
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coefficient vector, this is a time-independent effective Schrodinger equation 

He = Ec (4) 

where the effective renormahsed superblock Hamiltonian H satisfies E = = c''"Hc. 

The renormahsation tensors at each position are defined from the coefficient tensor at the 
same position, i.e. C" defines L" and R", via intermediate left and right density matrices. 
To obtain the left density matrix D^, we view the tensor C" as a matrix C indexed by 
{ln),r, where / is the row index of C", then = CC^ The right density matrix D/j is 
defined in a similar way, we view the tensor as a matrix C indexed by /, {nr), where r is 
the column index of C", and D/j = C^C. The renormahsation tensors L", R", when viewed 
as matrices L, R in the appropriate way, are obtained from the M eigenvectors (with largest 
weights) of the the density matrix Y)l and Dr respectively i.e. 

Dz,L = L(cri . . . CrM)diag, 0"! > . . . > (Tm, (L(/„)i = L"-) (5) 
DijR = R(ai . . . aA/)diag, (^l > <^2 ■ ■ ■ > CTM, (R-(rn)i = R^r) (6) 

More explicitly, writing the eigenvectors of the left and right density matrices as P, r*, 

Di.r = lV„ (71 > (72 > (73... (7) 
Brt' = rV„ CTi > (72 > • • • (8) 

L", R" are constructed by assigning the elements of the eigenvectors to the tensors in the 
following way 

LJ, = r(„,.), z = l...M (9) 
I^r, = r;„,), t = l...M (10) 



In Ref. |29|], we showed that satisfying the solution conditions Eqs. (jll), ([5]), ([6]) for 
C", L", R" is formally equivalent to minimising the DMRG energy subject to normalisation 
and the orthogonality constraints ([2]), ([3]). We can formally extend the DMRG theory to 
ime-dependent scenarios by making stationary the Dirac-Frenkel action {^{id/dt — Hl'^) 
12] subject to the same normalisation and orthogonality constraints. (Interestingly, the 
Dirac-Frenkel action has recently been independently rederived in the DMRG context in 



Ref. 30|). For the coefficient vector the time-evolution is then given by an effective time- 
dependent Schrodinger equation 

idtc = He (11) 
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The corresponding L" and R" remain defined by Eqs. ([5]), ([6]). 

III. COUPLED-PERTURBED DENSITY MATRIX RENORMALIZATION 

GROUP RESPONSE EQUATIONS 

We now consider the effect of an external perturbation. We start with a time-independent 
perturbation V. In the superblock basis {l}^{np}^{r}, this yields the effective perturbation 
V which satisfies (\1>|V|\P) = c+Vc. 

In response to this perturbation, the L", C", R" tensors each can be expanded in orders 
of \V\, giving 





1 + L"!^l + . . . 


(12) 




'] + c"[ii + . . . 


(13) 


R" = 


'1 + R"W + . . . 


(14) 



Thus the first-order DMRG wavefunction for the block-configuration in Fig. [T] takes the 
general form 

|v]/W) = ^ [(L"^[^l . . . C*''"' . . . R"*^) + . . . + (L"if°' . . . 0""^^^ . . . R"''!"!) 

{n} 

+ ... + (L"i[°l . . . C"-[°l . . . R"'=W)] \n^n2 ...rik) (15) 

We now derive the response equations satisfied by each of the quantities L"l^], C"[^],R"f^l 
These are obtained by the perturbation expansion of the solution conditions (jll), ([5]), IQ. 
For the coefficient vector, this yields 

(hM + AHW + VW + . . .)(c[°l + cW + ...) = (eM + + . . .)(cM + cW + . . .) (16) 

Note the first-order change in the Hamiltonian AH'^'. This arises because the effective 
Hamiltonian in the superblock basis H depends on the renormalisation tensors L", R" (which 
define the renormalised basis) and so first-order changes in those tensors lead to a first- 
order change in the effective Hamiltonian. (The construction of AH^^^ is described later in 
Sec. II V|) . Gathering first-order terms and enforcing intermediate normalisation through the 
projector Q = 1 — c^c'^l^ gives 

(hM - = -Q(AHW + vW)cM (17) 
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Because AH' ' depends on the first-order wavefunction through its dependence on the L", R" 
tensors, Eq. ( |T71) must be solved self-consistently. It is therefore a coupled-perturbed re- 
sponse equation, analogous to the coupled-perturbed orbital equations that arise in the 
Hartree-Fock theory of response. 

The first-order coefficients 0"^^^ define first-order renormalisation tensors at the same 
site L"W,R"W. Viewing C"'°1,C"W as a matrices in the appropriate fashion, we obtain 
first-order left and right density matrices 

dW = CMcWt + cWc[°lt, (C(„,),, = d) (18) 
dW = CM^CW + CWtclo], (c,(,,) = d) (19) 

In response to the change in the density matrices, the eigenvectors have a perturbation 
expansion 

r = + pi^i + . . . (20) 

= + r*[il + . . . (21) 

and we can set up corresponding response equations 

(d{;1 - a4)rW = -QlDWfM (22) 
(D^l - aa)r^[il = -Q^DWr^M (23) 

where Ql, Q_r project out the span of Di, T>R respectively, i.e. = Qi? = 

1 — r^l'^lr't"]"''. We assign the elements of each of the M perturbed vectors pl^^r^W 
according to Eq. ([9]), (ITO!) . to define L"l^], R"l^l The response equations for a time- dependent 
perturbation may be obtained in an analogous way as above. We consider for simplicity a 
perturbation with a single Fourier component, 

V{t) = l/e*^* + V*e-'^^ (24) 

We expand the L", C", R"" tensors in terms of orders of \V\, 

L"(t) = (L"M + L"W(t) + . . .)e-'^™* (25) 
C"(t) = (C"[°l + C"[^l (t) + . . .)e-'^'°'* (26) 
R"(t) = (R"M + R"[^l {t) + .. .)e-'^™* (27) 
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For the coefficient vector, we substitute this expansion into the effective time-dependent 
Schrodinger equation ffTTl) and identify terms with frequencies uj, —u, giving 

(HM - (^M + a;)l)cW(cu) = -Q(AHW(cu) + vW)cM (28) 
(hM - (E™ - cu)l)cW(-cu) = -Q(AHW(-cu) + vW*)c[°l (29) 

where Q is the projector defined in Eq. f|T7j) . The ffist-order frequency perturbed wavefunc- 
tions then define ffist-order perturbed density matrices Di(c<j), Di(—co'), 0/2(0;), D/j(—ct;), 
which can be used to obtain L"W(cj), L"W(-tu), R"W(tu), R"W(-cu) through Eqs. (|22|). (|23D. 

A. Response properties 

Once we obtain the ffist-order response of the DMRG wavefunction we can evaluate 
response properties of interest. We take as our example here the dipole-dipole response 
function or polarisability. For a uniform static electric field £i, the dipole moment is ex- 
panded as 

i^i = /^f ' + XI "^^'^J' + ■ ■ ■ ' hj ■■■ ^x,y,z (30) 
j 

which defines the static polarisability atj as the ffist-order change in the dipole moment. 
Within the DMRG response theory, the polarisability is therefore obtained as 

a,, = cMt;,f]cW + cWVfcM + c[°lt;,W.^c[°l (31) 

Here /ij is the effective dipole operator in the superblock basis, and c^^^ is the ffist-order 
wavefunction in response to an electric field in the j direction. Note the additional contri- 
bution jJ^jy This is the change in the effective dipole operator /ij due to the response of 
the L", R" tensors to an applied field in the j direction. This quantity is constructed in a 
similar way to the effective Hamiltonian AH'^'. 

For a frequency dependent electric field Si{t), we expand the dipole moment as 

fiiit) = /if due-'^'aij{u)£j{uj) + ... t,j...ex,y,z (32) 

j 

where o.ij{ijj) and £j{uj) are the u frequency components of the frequency dependent polar- 
isability and electric field. aij{uj) contains two contributions, one from the e*'^* component 
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of the applied perturbation, one from the e component. The final expression for aij{uj) 
therefore reads as 

aij{uj) = Gijiuj) + Gij{-uj) (33) 
G,,{i.) = cmf4>^cf{u) + cf (^)/.flcM + c[°VS)(^)c™ (34) 

Gij{u) and Gij{uj) are obtained from two separate response calculations, solving Eq. ( l28l) . 
(l29l) respectively. 



B. Comparison to other DMRG response theories 



So far we have derived a DMRG theory of response that was based on expanding the 
solution conditions satisfied by the DMRG wavefunction in terms of the applied perturba- 
tion. This corresponds to an analytic theory of response in the following way. Consider a 
time-independent perturbation for simplicity. Let us consider minimising the energy of the 
DMRG wavefunction, for some fixed number of states M , with respect to the full Hamilto- 
nian (with the perturbation) H = if '^l + W^^^ where A is used to scale the strength of the 
perturbation. This gives a wavefunction ^(A) and an energy E{X). The first-order wave- 
function and corresponding first-, second-, and third-order energies calculated with the 
analytic DMRG response theory correspond exactly to the following derivatives 

9^(A) 



dX 



A=0 



i^[3] 



dX 
d^E{X) 



A=0 



dX^ 
d^E{X) 



A=0 



A=0 



(35) 
(36) 
(37) 
(38) 



Analogous statements for time-dependent perturbations can be made by considering an 



appropriate quasi-energy 



HQ, 



19|. 



The analytic approach to DMRG response does not represent the only way to obtain 
response within the DMRG. Existing DMRG response methods use various related adaptive 
basis approaches, commonly in two categories, the Lanczos vector method and the 
dynamical density matrix renormalisation group 23j. The dynamical density matrix renor- 



malisation group is established as the most accurate approach to response properties and 



9 



we shall focus on it here. (Note the dynamical density matrix renormalisation group and 



correction vector methods 



21 



22 



25l | are essentially the same but differ in the algorithm 
used to solve the response equations. In fact, if the response quantities are evaluated using 
a quadratic functional of the correction vector such as Eq. ( l48i) . it is possible to obtain 
quadratic errors with the correction vector method without the explicit minimisation as 
used in the dynamical DMRG). 

In the dynamical DMRG the ansatz for the zeroth and first-order wavefunction are both 
modified relative to the unperturbed DMRG wavefunction, i.e. 

{n} 

\^m)=J2 L"^[°l . . . C"-W . . . R"^M \n1n2 ...Uk) (40) 

The tildes indicate that the L", C", R" tensors appearing in Eqs. (l39i) . even for the zeroth 
order wavefunction, do not correspond to the same tensors obtained in a DMRG calculation 
without the perturbation. The zeroth and first-order coefficient vectors are obtained from 
the effective Schrodinger equation (jlj) and an uncoupled response equation, e.g. 

(hM - (eM + o;)l)cW(^) = -QVWcM (41) 

The dynamical DMRG ansatz is able to capture the response of the L" and R" tensors in 
an average way, because it uses L", R" that are different from those in the unperturbed 
DMRG calculation. Specifically, the left and right renormalisation tensors at each block 
configuration are obtained as eigenvectors of modified left and right density matrices, where 
the density matrices corresponding to 0^,0^^^, v = V'^'c'^l are all averaged together i.e. for 

= aC^CMt + /5cWcWt + ^(Vc[°l)(VcM)t (42) 

where a + /5 + 7 = 1 and in the last term we are interpreting the perturbation multiplied 
by the zeroth order wavefunction V^^^cM as a matrix in the same way as c'^l is interpreted 
as a matrix. (Note that the above is for real frequencies; when for complex frequencies, 



one typically separates the imaginary and real contributions of the response vector 23|). 
Because the density matrix contains information on the perturbation and the response, 
the DMRG basis is "adapted" to the perturbation being considered. While this is very 
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simple to implement within a standard DMRG algorithm and has proven very successful, 
one drawback relative to the analytic response approach is that a single set of DMRG basis 
states is being used to represent several quantities, including both the zeroth order and 
response vectors. For this reason, we can expect some loss of accuracy with this method for 
small M calculations relative to the analytic response method. 



IV. IMPLEMENTATION 



We have implemented the analytic DMRG response theory as described above. This 
consists of three parts: solving the coupled-perturbed equation f|T7|) for the first-order co- 
efficient vector solving for the first-order renormalisation tensors L"^^], R"[^l ([5]), ([6]), 
and constructing the first-order effective Hamiltonian AH'^' and necessary intermediates, as 
well as other first-order operators needed for properties (e.g. in Eq. [3T1) . The first two 
parts are quite straightforward: we solve the coupled-perturbed equation ( |T7I) using a Krylov 
subspace iterative solver with preconditioning, and to obtain the first-order renormalisation 
tensors (122|) , (123!) we use explicit Rayleigh-Schrodinger expressions for the first-order density 
matrix eigenvectors 

p[i] = _[ ^ p[o] /43X 

j=M+l 

T-iMt-n[i]„i[o] 



j=M+l 

We now focus on the implementation to obtain AH^^' and related quantities such as 
We recall that the effective Hamiltonian H'"' is expressed as a tensor product of operators 
on the left and right blocks (we consider the single-site • in the block configuration Fig. [1] 
to be part of the left block for simplicity) 

H = ^^,,01®0^^ (45) 

where acts only the left block and Or acts only on the right block, and we assume that 
(8> takes into account the appropriate parity factors associated with the fermion character 
of the operators (see e.g. Ref. [2, 14])- The first-order Hamiltonian is constructed from the 
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response of the operators O^, Or, through 



AHW = 5^«;,,(Of of + of 'of ) 



(46) 



We therefore need to calculate the first-order operators O^ , O^ . These are built up sequen- 
tially through the blocking steps in the sweep much like the zeroth order operators. The 
renormalisation transformation TZ of the first-order operator at a given block configuration 
in a left— s>right sweep, is given by 



where we have used the underline to indicate that the operators refer to blocked operators 
(i.e. for the left block plus the single-site), and the renormalisation tensors are interpreted 
as matrices L as described in Eq. ([5]). At the beginning of the left— >right sweep, O^^ = for 
all such operators. Analogous expressions hold for the right— i>left sweep and the operators 



The full sweep algorithm for the DMRG analytic response can be summarised as follows: 

1. Converge a standard DMRG algorithm for the state of interest and store all interme- 
diate zeroth-order operators O^^, O^' and tensors L"'''!, C"'^^\ R"''''. 

2. Set all 0™,0™ = 

3. Start a sweep (left— >right) 

• Set all oK' to 

• At each block configuration: 

• Solve coupled perturbed response equation, Eq. ( ITTl) . AH^^' is constructed using 
current best guesses for O^^, O^' 

• Solve for perturbed density matrix eigenvectors and L"[^l, Eq. (1221) 

• Update all Oj^ using Eq. ^ 

4. Start a sweep (right— s>left), analogous to (left— i>right) sweep 

5. Loop to 3. until convergence. 



(47) 



O^. 
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FIG. 2: Oligo-di-acetylenes, with the long-axis moment of inertia ahgned with the x-coordinate. 
This is the axis along which the polarisabilities are evaluated. 

6. Evaluate response properties (e.g. as in Sec. IIII A\} 

We note that the cost of a single sweep for the analytic response has the same order of 
computational and storage cost as an ordinary sweep in the DMRG calculation, which, 
for the ab-initio Hamiltonian is 0{M^k^) +0(M^fc^) computation, 0{M'^k'^) memory, and 
0{M'^k^) disk, where k is the number of correlated orbitals. The memory cost is roughly 
twice that for the calculation of the energy because of storage of the first-order operators as 
well as the zeroth-order operators. 




V. STATIC AND FREQUENCY-DEPENDENT POLARIZABILITIES OF 

OLIGO-DI-ACETYLENES 



As an initial test of the analytic DMRG response theory and implementation, we have 
calculated static and frequency-dependent longitudinal polarisabilities of several oligo-di- 
acetylenes using the analytic DMRG response theory, the dynamical DMRG method, and 
the linear- response coupled cluster method. Long oligo-di-acetylenes are of interest due to 
their large third-order non-linear polarisability 31|- While we will calculate only the linear 
polarisability here, the same analytic derivative techniques can in principle be extended to 
higher order polarisabilities and non-linear optical response. 

We carried out calculations on short all-trans oligo-di-acetylenes (ODAs), 2-ODA CgHg, 4- 
ODA CigHio, 6-ODA C24H14. Optimised geometries were obtained at the density functional 
;heory B3LYP 32], 33] level in a correlation consistent Dunning double-zeta (cc-pVDZ) basis 
34| . Subsequent Hartree-Fock, DMRG, and coupled cluster (CC) calculations were carried 



out in a minimal ST0-6G Gaussian basis 



34 



35] . We realise that this basis is too small for 



the quantitative calculation of polarisabilities, but it has been chosen to enable a preliminary 
study. Also, we note that qualitative trends in polarisabihties can be captured using rather 
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small basis sets of split- valence quality [3l|]. The Hartree-Fock calculations were used to 
determine molecular orbitals with a and vr character. All a orbitals were kept frozen in the 
DMRG response calculations, and the vr orbitals were localised. Calculated polarisabilities 
refer to the Uxx component, where the x-axis is aligned with the long moment of inertia axis 
of the molecules (see Fig. [2]). The DMRG response calculations used an active space oi pz 
orbitals only, corresponding to an (8e, Sorb) active space for 2-ODA, a (16e, 16orb) active 
space for 4-ODA, and a (24e, 24orb) active space for 6-ODA. For the analytic response 
DMRG calculations using M states we first converged a ground-state DMRG calculation 
with M states using the one-site algorithm, and used this as the starting point for the 
response calculation. 

In addition to the analytic response DMRG calculations, we carried out calculations 
using the dynamical DMRG method for comparison. The dynamical DMRG polarisabilities 
were obtained by solving the linear response equation in the dynamical DMRG basis {uol — 
jj[o])c[i]). = Q^jC^^] just as in the correction vector method, but the resulting polarisabilities 
were evaluated using the quadratic functional 

G., = c^^\ul - H[°l)cW + cMt^^cW + cf^x^cM J (48) 



which ensures that the obtained polarisability is quadratic in the error in ct^l [361, ISTj , which is 
the hallmark of the dynamical DMRG approach. For comparison, we also computed linear- 
response restricted coupled cluster polarisabilities at the singles and doubles level [l^ , both 
at the all electron level, and within the active space only, using the Psi3 [38] package. 

We note one issue that arises with the response DMRG calculations in our initial imple- 
mentation as opposed to ordinary ground-state DMRG calculations. In ground-state DMRG 
calculations with the one-site algorithm, we are generally able to converge the DMRG en- 
ergy from sweep to sweep to very high accuracy, e.g. nanoHartrees. However, in our initial 
response implementation, we were not able to converge the calculated polarisabilities to 
similar accuracy. Typically the forward and backwards sweeps would converge to somewhat 
different results, and even between consecutive forwards (or backwards) sweeps, the polar- 
isability would oscillate somewhat. This was true both for the dynamical DMRG and the 
analytic response DMRG calculations. The oscillation can be quite severe, particularly for 
small M calculations and for higher frequencies that are nearer to a pole (e.g. at frequency 
ui = 0.2 a.u.) and reflects the greater sensitivity of the response calculation to the discarded 
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states in the density matrix. In our results, we report the average polarisabihty of the last 
4 sweeps, together with twice the standard deviation. These results are reported in tabled 

From table [T] we make the following observations about the relative performance of the 
analytic DMRG response method relative to the dynamical DMRG method that has been 
commonly used. For small M (e.g. M=25) the analytic DMRG response method is clearly 
superior. Whereas the dynamical DMRG method produces poor polarisabilities for M=25, 
in error by more than 50% in some cases, the analytic DMRG polarisabilities are quite 
reasonable at M=25 and typically in error by less than 1%. This is consistent with our 
discussion in section IIII Bl where we argue the the dynamical DMRG method suffers from 
using the same set of DMRG basis states to represent both the zeroth order DMRG vector 
as well as the response and perturbation vectors. Thus, for small M there simply are not 
enough DMRG states to yield a meaningful result in the dynamical DMRG. Both methods 
converge as M increases. For the most accurate calculations (M=250), although both 
methods perform well, the dynamical DMRG polarisabilities appear slightly better than 
the analytic DMRG polarisabilities. However, this appears to be related to the instabilities 
in the convergence of the analytic DMRG response sweeps; whereas the oscillations in the 
dynamical DMRG sweeps vanish for larger M, they still remain for the analytic DMRG 
sweeps. From the 2(T values, we see that currently we can only conclude that the analytic 
and dynamical DMRG response methods are comparable for larger M. 

Observing the trends in the polarisabilities, we see that the polarisabilities increase as 
the applied frequency increases which is what one would expect since we are approaching 
the first excitonic ^Bu pole. We are not able to converge our response calculations very close 
to a pole because of the large norm in The standard solution to this is to include a 
small imaginary broadening in uj. However, a straightforward incorporation of broadening 
leads to complex operators in the analytic theory which we have not yet implemented. 

It is often the case that one wishes to determine an entire spectrum, i.e. some response 
property for a very large range of uj. While in the dynamical DMRG this is usually per- 
formed by scanning through uj (with some small imaginary component) and performing a 
response calculation for each frequency, it may be more appropriate in the analytic response 
approach to adopt a different strategy. The coupled-perturbed response equations may be 
viewed as a linear eigenvalue problem for the excitation energies (i.e. poles) and may be 
solved in this way, in the same way that the time- dependent Hartree-Fock or time- dependent 
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TABLE I: Static and frequency dependent polarisabilities in a.u. of oligo-di-acetylcnes, with 2, 
4, 6 monomers (2-ODA, 4-ODA, 6-ODA). D stands for dynamical DMRG, A stands for analytic 
response theory, u is the frequency (in a.u.), M refers to the number of states in the DMRG 
wavefunction. The numbers in brackets do not represent intrinsic truncation error from 
finite M but represent the numerical convergence of the DMRG sweep, since the forwards and 
backwards sweeps typically converge to slightly different results. The bracketed number is twice 
the standard deviation {2a) in the last 4 forward and backwards sweeps. See text for further 
discussion. 





2-ODA 4-ODA 6-ODA 


w M 


D A 


D {2a) A {2a) 


D {2a) A {2a) 


0.00 25 

50 

250 

1000 

LR-CCSD 


52.77 52.89 
52.89 52.89 
52.88 52.88 
n.a. n.a. 
53.38 


144.16 (0.03) 145.21 (0.04) 
146.07 (0.01) 145.74 (0.09) 

145.75 (0.01) 145.80 (0.01) 
145.77 (0.01) 145.81 (0.00) 
148.15 


354.28 (17.96) 243.65 (0.06) 
246.04 (0.02) 245.06 (0.07) 

245.20 (0.00) 245.27 (0.03) 
245.13 (0.10) 245.14 (0.02) 
249.67 


0.05 25 
50 

250 

LR-CCSD 


53.98 53.96 
54.07 54.07 

54.06 54.07 
54.62 


148.46 (0.02) 149.80 (0.04) 
150.64 (0.01) 150.26 (0.07) 
150.37 (0.00) 150.39 (0.04) 
153.19 


449.82 (35.15) 252.00 (0.14) 
254.61 (0.02) 253.62 (0.13) 
253.87 (0.00) 253.92 (0.02) 
259.40 


0.10 25 
50 
250 

LR-CCSD 


57.83 57.57 
57.99 57.99 
57.99 58.00 

58.72 


163.62 (0.03) 165.42 (0.13) 
166.46 (0.02) 166.11 (0.05) 
166.19 (0.00) 166.23 (0.02) 
170.76 


462.00 (22.55) 282.05 (0.25) 
284.81 (0.03) 283.96 (0.22) 
284.30 (0.00) 284.26 (0.21) 
294.16 


0.15 25 
50 
250 

LR-CCSD 


65.85 64.97 
66.07 66.06 
66.05 66.08 
67.22 


195.14 (0.07) 201.06 (0.17) 
202.51 (0.03) 202.24 (0.09) 
202.45 (0.00) 202.49 (0.04) 
212.20 


557.18 (114.72) 353.66 (0.57) 
357.02 (0.05) 356.37 (0.20) 
357.26 (0.00) 357.10 (0.10) 
381.68 


0.20 25 
50 
250 

LR-CCSD 


82.03 79.89 
82.57 82.54 
82.56 82.60 
84.83 


279.06 (0.35) 294.06 (0.89) 
296.86 (0.62) 295.83 (1.67) 
296.71 (0.55) 296.44 (0.06) 
328.71 


520.61 (84.68) 564.50 (1.38) 
564.25 (16.84) 566.94 (0.89) 
571.44 (0.71) 571.63 (1.73) 
682.10 
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2 3 4 5 6 



number of monomers 



FIG. 3: Scaling of total and active space polarisabilities per monomer. 

density functional equations are solved as an eigenvalue problem to obtain excitation ener- 
gies. Once a sufficient number of poles are obtained, the spectrum can then be reconstructed 
analytically. 

Comparing the DMRG polarisabilities and the coupled cluster polarisabilities, we see that 
the coupled cluster polarisabilities are generally quite good even at the singles and doubles 
level. (They appear to consistently overestimate the polarisability by only a few percent). 
This is not surprising since by virtue of the one-electron nature of the dipole operator, 
the linear polarisability only samples states with single-excitation character relative to the 
ground-state. Such excited states are well captured by CCSD theory. However, earlier 
studies indicate that the overall spectrum in conjugated systems (including e.g. doubly 
excited and triplet excited states) is poorly reproduced by coupled cluster theory jsot, and 
so we would expect much larger discrepancies between the CC and DMRG description of 
third-order non-linear optical response. 

In Fig. [3] we plot the static active space and total polarisabilities {u = 0) per monomer 
calculated using the analytic DMRG response theory as a function of the number of di- 
acetylene monomers in the calculation. The total polarisability for the DMRG calculations 
is obtained using the core-correction from the linear-response coupled cluster calculations 
i.e. 

^ tot ^tot ^ act I ^ act fAn\ 

"DMRG = "cC - "cC + "DMRG (49) 

We see a slow saturation of the polarisability per monomer as a function of the chain length. 
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although the polarisabihty is not yet fully saturated at the 6-ODA level. While larger 
basis set calculations and calculations on longer chains are necessary to obtain a definitive 
conclusion, we note that our results are consistent with early semi-empirical calculations 
which indicate an onset of saturation between 2-ODA and 3-ODA {4^. 

VI. CONCLUSIONS 

In the current work we have described an analytic approach to the calculation of response 
quantities in the density matrix renormalisation group. The analytic response method is 
familiar from other electronic structure theories but has not so far been developed within 
the density matrix renormalization group. The analytic response implementation does not 
change the computational cost of the ground-state DMRG calculation by more than a con- 
stant factor. Compared to the popular dynamical density matrix renormalisation group 
approach we find that the analytic response method produces considerably more accurate 
response quantities when using a small number of DMRG states, without any greater com- 
putational cost. While it is simpler within the dynamical DMRG to implement higher-order 
response properties and complex frequencies, based on our investigations, the improved ac- 
curacy of the analytic response approach may justify the additional implementation effort. 
In future work, we will explore both higher-order response quantities and determination of 
complete spectra using the analytic DMRG response approach. 
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